Cluster growing process and a sequence of magic numbers. 
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We present a new theoretical framework for modelling the cluster growing process. 
Starting from the initial tetrahedral cluster configuration, adding new atoms to the 
system and absorbing its energy at each step, we find cluster growing paths up to 
the cluster sizes of more than 100 atoms. We demonstrate that in this way all known 
global minimum structures of the Lennard-Jonnes (LJ) clusters can be found. Our 
method provides an efficient tool for the calculation and analysis of atomic cluster 
structure. With its use we justify the magic numbers sequence for the clusters of 
noble gases atoms and compare it with experimental observations. We report the 
striking correspondence of the peaks in the dependence on cluster size of the second 
derivative of the binding energy per atom calculated for the chain of LJ-clusters 
based on the icosahedral symmetry with the peaks in the abundance mass spectra 
experimentally measured for the clusters of noble gases atoms. Our method serves an 
efficient alternative to the global optimization techniques based on the Monte-Carlo 
simulations and it can be applied for the solution of a broad variety of problems in 
which atomic cluster structure is important. 

It is well known that the sequence of cluster magic numbers carries essential information 
about the electronic and ionic structure of the cluster Understanding of the the clus- 
ter magic numbers is often equivalent or nearly equivalent to the understanding of cluster 
electronic and ionic structure. A good example of this kind is the observation of the magic 
numbers in the mass spectrum of sodium clusters 0. In this case, the magic numbers were 
explained by the delocalised electron shell closings (see || and references therein). Another 
example is the the discovery of fullerenes, and in particular the Cqq molecule |3J] , which was 
made by means of the carbon clusters mass spectroscopy. 
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The formation of a sequence of cluster magic numbers should be closely connected to 
the mechanisms of cluster formation and growing. It is natural to expect that one can 
explain the magic numbers sequence and find the most stable cluster isomers by modelling 
mechanisms of cluster assembling and growing. On the other hand, these mechanisms are 
of interest on their own and the correct sequence of the magic numbers found in such a 
simulation can be considered as a proof of validity of the cluster formation model. 

The problem of magic clusters is closely connected to the problem of searching for global 
minima on the cluster multidimentional potential energy surface. The number of local 
minima on the potential energy surface increases exponentially with the growth cluster size 
and is estimated to be of the order of 10 43 for N = 100 Q. Thus, searching for global minima 
becomes increasingly difficult problem for large clusters. There are different algorithms and 
methods of the global minimisation, which have been employed for the global minimisation 
of atomic cluster systems (see |]J and references therein). These techniques are often based 
on the Monte-Carlo simulations. 

The algorithm which we describe in this work is based on the dynamic searching for the 
most stable cluster isomers in the cluster growing process. Our calculations demonstrate 
that our approach is an efficient alternative to the known techniques of the cluster global 
minimisation. The big advantage of our approach consists in the fact that it allows to study 
not just the optimized cluster geometries, but also their formation mechanisms. 

In the present work we approach the formulated problem in a most simple, but general 
form. In our most simple scenario, we assume that atoms in a cluster are bound by Lennard- 
Jones potentials and the cluster growing takes place atom by atom. In this process, new 
atoms are placed on the cluster surface in the middle of the cluster faces. Then, all atoms 
in the system are allowed to move, while the energy of the system is decreased. The motion 
of the atoms is stopped, when the energy minimum is reached. The geometries and energies 
of all cluster isomers found in this way are stored and analysed. The most stable cluster 
configuration (cluster isomer) is then used as a starting configuration for the next step of 
the cluster growing process. 

Starting from the initial tetrahedral cluster configuration and using the strategy described 
above, we have analysed cluster growing paths up to the cluster sizes of more than 100 
atoms. We have found that in this way practically all known global minimum structures 
of the Lennard-Jonnes clusters (see |1[ and references therein) can be determined, which 
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proves that our method is indeed an efficient alternative to other cluster global optimization 
techniques such as basin hoping algorithm [|TJ. 

In our model we consider an atomic cluster as a group of atoms that interact with each 
other by pairing forces. The interaction potential between two atoms in the cluster can, in 
principle, be arbitrary. In this work, we use the Lennard- Jones (LJ) potential: 



where r is the interatomic distance, — e is the depth of the potential well (e > 0), 2 1 / 6 cx is 
the pair bonding length. 

The constants in the potentials allow one to model various types of clusters for which 
LJ paring force approximation is reasonable. The most natural systems of this kind are the 
clusters consisting of noble gases atoms Ne, Ar, Kr, Xe. The magic numbers for this type 
of clusters have been experimentally determined in [JSj, |J. In our modelling of the cluster 
growing process we focus on this example and consider it below in detail. The constants in 
the L J potential appropriate for the noble gases atoms one can find in . The L J forces are 
also appropriate for modelling nuclear clusters consisting of alpha particles M. Note that 
for the L J clusters it is always possible to chose the coordinate scale so that a = 1 . It makes 
all LJ cluster systems scalable. They differ only by the choice of the energy parameter e 
and the mass of a single constituent (atom). 

In our approach the atomic motion in the cluster is described by the Newton equations 
with the LJ pairing forces. The system of coupled equations for all atoms in the cluster 
are solved numerically using the 4-th order Runge-Kutta method. The primary goal in 
this simulation was to find the solutions of the equations that lead to the stable cluster 
configurations and then to chose energetically the most favourable one. The choice of initial 
conditions for the simulation and the algorithm for the solution of this problem are described 
below. 

Our cluster searching algorithm is constructed on the idea that each minimum on the 
cluster potential energy surface corresponds to the situation, when all the atoms are located 
in their equilibrium positions. A minimum can be found by allowing atoms to move, starting 
from a certain initial cluster configuration, and by absorbing all their kinetic energy in the 
most efficient way. If the starting cluster configuration for N + 1 atoms has been chosen on 
the basis of the global minimum structure for iV atoms, then it is natural to assume, and 
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we prove this in the present work, that the global minimum structure for N + 1 atoms can 
be easily found. The success of this procedure reflects the fact that in nature clusters in 
their global minima often emerge namely in the cluster growing process, which we simulate 
in such calculation. 

We have employed the following algorithm for the kinetic energy absorption. At each 
step of the calculation we consider the motion of one atom only, which undergoes the action 
of the maximum force. At the point, in which the kinetic energy of the selected atom is 
maximum, we set the absolute value of its velocity to zero. This point corresponds to the 
minimum of the potential well at which the selected atom moves. When the selected atom 
is brought to the equilibrium position, the next atom is selected to move and the procedure 
of the kinetic energy absorption repeats. The calculation stops when all the atoms are in 
equilibrium. 

We have considered a number of scenario of the cluster growing on the basis of the 
developed algorithm for finding the stable cluster configurations. 

In the most simple scenario clusters of N + 1 atoms are generated from the N-atomic 
clusters by adding one atom to the system. In this case the initial conditions for the sim- 
ulation of N+l-atomic clusters are obtained on the basis of the chosen N-atomic cluster 
configuration by calculating the coordinates of an extra atom added to the system on a cer- 
tain rule. We have probed the following paths: the new atom can be added either (Al) to 
the center of mass of the cluster, or (A2) randomly outside the cluster, but near its surface, 
or [A3) to the centers of mass of all the faces of the cluster (here, the cluster is considered 
as a polyhedron), or (A4) to the points that are close to the centers of all the faces of the 
cluster, located from both sides of the face on the perpendicular to it, {A5) to the centers 
of mass of the faces laying on the cluster surface. 

The choice of the method how to add atoms to the system depends on the problem to be 
solved. For example, with Al and A2 methods large clusters consisting of many particles can 
be generated rather quickly. The A2 method is especially fast, because adding one atom to 
the boundary of the cluster usually does not lead to the recalculation of its central part. The 
A3 axxi A4 methods can be used for searching the most stable, i.e. energetically favourable, 
cluster configurations or for finding cluster isomers with some other specific properties. The 
A4 method leads to finding more cluster isomers than the A3 one, but it takes more CPU 
time. The A5 method is especially convenient for modelling the cluster growing process 
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which we focus on in this paper. Using this method one can generate the cluster growing 
paths for the most stable cluster isomers. 

When considering the cluster growing process, new atoms should be added to the system 
starting from the initially chosen cluster configuration step by step until the desired cluster 
size is reached. Each new step of the cluster growing should be made with the use of the 
methods A1-A5. The criteria for the cluster selection in this process can be as follows: at 
every step (SE1) one of the clusters with the minimum number of atoms is selected, or 
(SE2) the cluster with the minimum energy among the already found stable clusters of the 
maximum size is selected, or (SE3) the cluster with the maximum energy among the already 
found stable clusters of the maximum size is selected. 

The SE1 criterion is relevant in the situation, when the full search of cluster isomers is 
needed. It is applicable to the systems with relatively small number of particles. The SE2 
criterion is relevant for modelling the cluster growing process. It turns out to be very efficient 
and leads to finding the most stable cluster configurations for a given number of particles. 
The SE3 criterion might be useful for the redirection of the cluster growing process towards 
the lower energy cluster isomers branches. 

Calculations performed with the use of the methods described above show that often 
clusters of higher symmetry group possess relatively low energy. Thus, the symmetric cluster 
configurations are often of particular interest. The process of searching the symmetric cluster 
configurations can be speed up significantly, if one performs the cluster growing process with 
the imposed symmetry constraints. This means that for obtaining a symmetric iV atomic 
cluster isomer from the initially chosen symmetric (N — M )-atomic configuration one should 
add M atoms to the surface of this isomer symmetrically. 

Using our algorithms we have examined various paths of the cluster growing process and 
determined the most stable isomers up to the cluster sizes of more than 100 atoms. The 
binding energies per atom as a function of cluster size for the calculated cluster chains are 
shown in figure [l|. In the insertion to figure |I| we present the experimentally measured 
abundance mass spectrum for the Ar clusters at 77K |5|. |(| . 

We have generated the chains of clusters based on the icosahedral, octahedral, tetrahedral 
and decahedral symmetries with the use of the A3-A5 and SE2-SE3 methods. In a few 
particular cases for N > 70, we have also used manual modifications of the starting cluster 
geometries. In all our calculations we have used the dimensionless form of the LJ-potential, 
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FIG. 1: Binding energy per atom for LJ-clusters as a function of cluster size calculated for the 
cluster chains based on the icosahedral, octahedral, tetrahedral and decahedral symmetry. In the 
insertion we present the experimentally measured abundance mass spectrum for the Ar clusters at 
77K U 

i.e. put a — 1. The potential constant has been chosen as e = 1/4. Such a choice of 
constants is the most universal. It allows one to rescale easily all the results to any concrete 
choice of a and e. 

Figure [i] shows that the most stable clusters are obtained on the basis of the icosahedral 
symmetry configurations with exceptions for N = 38, 75 < N < 77 and N = 98. In these 
cases the octahedral cluster symmetry becomes more favourable. The cluster chains based 
on the tetrahedral and decahedral symmetries have no intersections with the icosahedral 
chain of clusters although there is interplay between these two curves and the octahedral 
one. 

The main trend of the energy curves plotted in figure [l] can be understood on the basis 
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of the liquid drop model, according to which the cluster energy is the sum of the volume 
and the surface energy contributions: 

E N = -X V N + X S N 2/3 - XrN 1 / 3 (2) 

Here the first and the second terms describe the volume, and the surface cluster energy 
correspondingly The third term is the cluster energy arising due to the curvature of the 
cluster surface. Choosing constants in @ as Ay = 0.71703, As = 1.29302 and \r = 0.56757, 
one can fit the global energy minimum curve plotted in figure [l] with the accuracy less 
than one per cent. The deviations of the energy curves calculated for various chains of 
cluster isomers from the liquid drop model (§) are plotted in figure |2|. The curves for the 
icosahedral and the global energy minimum cluster chains go very close with each other and 
the peaks on these dependences indicate the increased stability of the corresponding magic 
clusters. The ratio between the volume and surface energies in (0) can be characterised by 
the dimensionless parameter 5 = Ay/As, being equal in our case to 5 = 0.5545. 

The dependence of the binding energies per atom for the most stable cluster configu- 
rations on N allows one to generate the sequence of the cluster magic numbers. In the 
insertion to figure |2| we plot the second derivatives AE 2 for the chain of icosahedral iso- 
mers. We compare the obtained dependence with the experimentally measured abundance 
mass spectrum for the Ar clusters at 77K ^|, (see insertion to figure [I]) and establish 
the striking correspondence of the peaks in the measured mass spectrum with those in the 
AE 2 dependence. Indeed, the magic numbers determined from A 2 E^ are in a very good 
agreement with the numbers experimentally measured for the Ar and Xe clusters: 13, 19, 
23, 26, 29, 32, 34, 43, 46, 49, 55, 61, 64, 71, 74, 81, 87, 91, 101, 109, 116, 119, 124, 131, 
136, 147 ||, ||. The most prominent peaks in this sequence 13, 55 and 147 correspond to 
the closed icosahedral shells, while other numbers correspond to the filling of various parts 
of the isosahedral shell. 

The connection between the second derivatives A 2 E^ and the peaks in the abundance 
mass spectrum of clusters one can understand using the following simple model. Let us 
assume that the mass spectrum of clusters is formed in the evaporation process. This 
means that changing the number of clusters, n^, of the size N in the cluster ensemble 
takes place due to the evaporation of an atom by the clusters of the size N and N + 1, i.e. 
An N ~ n N+ iWN+i~>N — u nWn^n-i, where the evaporation probabilities are proportional to 



8 




FIG. 2: Energy curves deviations from the liquid drop model (|2|) calculated for various cluster 
isomers chains. In the insertion we plot the second derivative A 2 En = E^+i + Ejq-i — 2En 
calculated for the icosahedral cluster isomers chain. 

E N +E 1 -E N+l E N _ 1 +E 1 -E N 

Wn+i^n ~ e kT and Wn-*n-i ~ e . Here T is the cluster temperature, 

k is the Bolzmann constant. In the limit AE^/kT <C 1, one derives An^ ~ ^aK-^Wi + 
Ejq^\ — 2E^jkT ~ A 2 En- These estimates demonstrate that the positive second derivative 
A 2 En should lead to the enhanced abundance of the corresponding clusters. 

In figure ^ we plot images of the magic clusters up to iV = 71. For iV = 32 and iV = 34, 
we present the icosahedral isomer and the one possessing the global energy minimum. We 
also plot the image of the octahedral iV = 38 cluster, which is found to be more stable than 
the clusters from the icosahedral chain. Experimentally iV = 38 is not found to be the magic 
cluster, although it is the global minimum cluster, being magic for the octahedral cluster 
chain (see figures || and §). This fact can be understood if one takes into account that 
different symmetry cluster chains are formed independently and the transition of clusters 
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FIG. 3: Geometries of the magic LJ-clusters. The labels indicate the cluster size and the cluster 
point symmetry group. 

from one chain to another at certain N is not possible. It is clear from the binding energy 
analysis that the icosahedral chain of clusters should be dominating. In experiments, clusters 
of the icosahedral chain mask clusters of other symmetry chains even in the situations when 
these other clusters are energetically more favourable, like it occurs for N = 38. 

In this paper we have discussed the classical models for the cluster growing process, but 
our ideas can be easily generalized on the quantum case and be applied to the cluster systems 
with different than LJ type of the inter-atomic interaction. It would be interesting to see 
to which extent the parameters of inter-atomic interaction can influence the cluster growing 
process and the corresponding sequence of magic numbers or whether the crystallization 
in the nuclear matter consisting of alpha particles and/or nucleons is possible. Studying 
cluster thermodynamic characteristics with the use of the developed technique is another 
interesting problem which is left opened for future considerations. 

The authors acknowledge support from the the Alexander von Humboldt Foundation and 
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